!> Solving the ODE y' = y - t*y with y(0) = 1 using implicit forward Euler
program implicit_forward_euler

    use display_module, only: display
    implicit none

    integer, parameter :: n = 10   ! number of time steps
    integer :: i, j, k             ! loop indices
    real :: t(n), y(n), f(n), dy(n), eps, h

    ! set initial values
    t(1) = 0.0
    y(1) = 1.0
    eps = sqrt(epsilon(1.0))
    h = 0.1

    do i = 2, n
        ! time step
        t(i) = t(i - 1) + h

        ! initial guess for y(i)
        y(i) = y(i - 1)

        ! not set the maximum number of iterations this time
        do
            ! calculate f(y(i)), y(i) = y(i-1) + h*f'(t(i), y(i))
            f(i) = y(i) - y(i - 1) - h*func(t(i), y(i))  ! implicit forward euler

            ! calculate dy/dy(i)
            dy(i) = 1.0 - h*dfunc_dy(t(i), y(i))  ! implicit forward euler

            ! check convergence
            if (abs(f(i)) < eps) exit

            ! update y(i)
            y(i) = y(i) - f(i)/dy(i)  ! newton iteration
        end do
    end do

    call display(t, 't:', .false.)
    call display(y, 'y:', .false.)

contains

    ! define the function
    function func(t, y) result(f)
        real, intent(in) :: t, y
        real :: f

        f = y - t*y

    end function func

    ! define the derivative of the function with respect to y
    function dfunc_dy(t, y) result(df_dy)
        real, intent(in) :: t, y
        real :: df_dy

        df_dy = 1.0 - t

    end function dfunc_dy

end program implicit_forward_euler
!> [vector: 10] t:
!>  0.000E+00,  1.000E-01,  2.000E-01,  3.000E-01,  4.000E-01,  5.000E-01,  6.000E-01,  7.000E-01,  8.000E-01,  9.000E-01
!> [vector: 10] y:
!>  1.000E+00,  1.099E+00,  1.194E+00,  1.284E+00,  1.366E+00,  1.438E+00,  1.498E+00,  1.545E+00,  1.576E+00,  1.592E+00
